STAD70

Q1

Suppose the four variables being used are denoted by \(X_1, \ldots, X_4\). Use the sample covariance matrix to estimate the variance of \(0.5 X_1+0.3 X_2+\) \(0.2 X_3\).

COV = var(Berndt)    # variance-covariance matrix of returns
w = c(.5, .3, .2, 0) # vector of weights

t(w) %*% COV %*% w   # variance of linear combination of returns (portfolio)

##             [,1]
## [1,] 0.004408865

Fit a multivariate- \(t\) model to the data using the function cov.trob in the MASS package. This function computes the MLE of the mean and covariance matrix with a fixed value of \(\nu\). To find the MLE of \(\nu\), the following code computes the profile log-likelihood for \(\nu\).

library(MASS) # needed for cov.trob
library(mnormt) # needed for dmt
df = seq(2.5, 8, 0.01)
n = length(df)
loglik_profile = rep(0, n)
for (i in 1:n) {
  fit = cov.trob(Berndt, nu = df[i])
  mu = as.vector(fit$center)
  sigma = matrix(fit$cov, nrow = 4)
  loglik_profile[i] = sum(log(dmt(Berndt, mean = fit$center, # mt for multivar t
    S = fit$cov, df = df[i])))
}
plot(loglik_profile, type = 'l')

Q2

Using the results produced by the code above, find the MLE of \(\nu\) and a 90% profile likelihood confidence interval for \(\nu\). Include your R code with your work. Also, plot the profile log-likelihood and indicate the MLE and the confidence interval on the plot.

Section 7.13.3 demonstrates how the MLE for a multivariate t-model can be fit directly with the optim function, rather than profile likelihood.

plot(df, 2 * loglik_profile, type = "l", cex.axis = 1.5,  # plot 2*log(likelihood)
     cex.lab = 1.5, ylab = "2 * log(likelihood)", lwd = 2)
abline(h = 2 * max(loglik_profile) )

# plot maximum
max_ind = which.max(loglik_profile) # Find index of maximum log(likelihood)
points( df[max_ind], 2*loglik_profile[max_ind], pch=15, cex=2)  
text( df[max_ind], 2*loglik_profile[max_ind], pos=1, labels = "MLE") 

# plot 90% CI level 
abline(h = 2 * max(loglik_profile) - qchisq(0.9, 1)) 
# find indexes of lower & upper bounds of 90% CI
lower_ind = head( which( loglik_profile > max(loglik_profile) - qchisq(0.9, 1 )/2 ), 1 ) 
upper_ind = tail( which( loglik_profile > max(loglik_profile) - qchisq(0.9, 1 )/2 ), 1 ) 
abline(v = df[lower_ind] )
abline(v = df[upper_ind] )
points( df[lower_ind], 2*loglik_profile[lower_ind], pch=16, cex=2)  # plot lower bound
text( df[lower_ind], 2*loglik_profile[lower_ind], pos=1, labels = "lower 90% \n CI bound") 
points( df[upper_ind], 2*loglik_profile[upper_ind], pch=16, cex=2)  # plot upper bound
text( df[upper_ind], 2*loglik_profile[upper_ind], pos=1, labels = "upper 90% \n CI bound") 

print( rbind( c( "lower", "max", "upper"), df[ c(lower_ind, max_ind, upper_ind)] ))

##      [,1]    [,2]   [,3]   
## [1,] "lower" "max"  "upper"
## [2,] "3.23"  "4.55" "6.59"

Q3

The following code generates and plots four bivariate samples. Each sample has univariate marginals that are standard t3-distributions. However, the dependencies are different.

library(MASS) # need for mvrnorm
par(mfrow=c(2,2))
N = 2500
nu = 3

set.seed(5640)
cov=matrix(c(1, 0.8, 0.8, 1), nrow = 2) # rho = 0.8 (correlated)
x= mvrnorm(N, mu = c(0, 0), Sigma = cov)
w = sqrt(nu / rchisq(N, df = nu))
x = x * cbind(w, w) # scale by same factor (tail dep)
plot(x, main = "(a)")

set.seed(5640)
cov=matrix(c(1, 0.8, 0.8, 1),nrow = 2) # rho = 0.8
x= mvrnorm(N, mu = c(0, 0), Sigma = cov)
w1 = sqrt(nu / rchisq(N, df = nu))
w2 = sqrt(nu / rchisq(N, df = nu))
x = x * cbind(w1, w2) # scale by dif factors
plot(x, main = "(b)")

set.seed(5640)
cov=matrix(c(1, 0, 0, 1), nrow = 2) # rho = 0
x= mvrnorm(N, mu = c(0, 0), Sigma = cov)
w1 = sqrt(nu / rchisq(N, df = nu))
w2 = sqrt(nu / rchisq(N, df = nu))
x = x * cbind(w1, w2) # scale by dif factors
plot(x, main = "(c)")

set.seed(5640)
cov=matrix(c(1, 0, 0, 1), nrow = 2) # rho = 0
x= mvrnorm(N, mu = c(0, 0), Sigma = cov)
w = sqrt(nu / rchisq(N, df = nu))
x = x * cbind(w, w) # scale by same factor
plot(x, main = "(d)")

Which sample has independent variates? Explain your answer.

• Samples (a) and (b) are (positively) correlated, because they both arise from (randomly scaled) correlated Normals (with \(\rho=0.8\)).
• Sample (d) exhibits tail-dependence, because the uncorrelated/independent Normals are scaled by the same random (\(\chi^2\)-based) factor.
• Only sample (c) is independent, because the two independent Normal variates are scaled by independent random factors.

Q4

Which sample has variates that are correlated but do not have tail dependence? Explain your answer.

• Sample (b) has correlated variates (because it is based on correlated Normals), but does not have tail-dependence, because its variates were scaled by independent random factors. The plot shows that extreme values seem to occur independently in the two dimensions (i.e., extreme values are either horizontally or vertically extreme, but not both.)

Q5

Which sample has variates that are uncorrelated but with tail dependence? Explain your answer.

• Sample (d) has uncorrelated variates (based on uncorrelated Normals) but with tail-dependence, because they are scaled by the same random factor. The plot shows that extreme values tend to occur simultaneously in both dimensions.

Q1

Consider the R code above.

What type of copula model has been sampled? Give the copula family, the correlation matrix, and any other parameters that specify the copula.

It is a t copula. The correlation matrix is “unstructured” meaning that it is an arbitrary correlation matrix. In fact, the correlation matrix is

\[ \left(\begin{array}{ccc} 1 & -0.6 & 0.75 \\ -0.6 & 1 & 0 \\ 0.75 & 0 & 1 \end{array}\right) \] The degrees of freedom parameter is 1.

Q2

Examine the scatterplot matrix (generated by line 6) and answer the questions below. Include the scatterplot matrix with your answer.

1. Components 2 and 3 are uncorrelated. Do they appear independent? Why or why not?

   Components 2 and 3 would be uniformly scattered over the unit square if they were independent.

   Clearly, the scatter is not uniform, so they do not appear independent.

2. Do you see signs of tail dependence? If so, where?

   The non-uniformity mentioned in (a) is that there are more data in the corners, which shows that extreme values tend to occur together, although because of the zero correlation, a positive extreme value of one component is equally likely to be paired with a positive or negative extreme value of the other component.

3. What are the effects of dependence upon the plots?

   The effects of tail dependence is the tendency of extreme values to pair. The negative correlation of components 1 and 2 shows in the concentration of the data along the diagonal from upper left to lower right. Positive extreme values in one component tend to pair with negative extreme values of ther other component.

   The positive correlation of components 2 and 3 shows in the concentration of the data along the diagonal from lower left to upper right. Positive extreme values in one component tend to pair with positive extreme values of ther other component.

4. The nonzero correlations in the copula do not have the same values as the corresponding sample correlations. Do you think this is just due to random variation or is something else going on? If there is another cause besides random variation, what might that be? To help answer this question, you can get confidence intervals for the Pearson correlation: For example, cor.test(rand_t_cop[,1],rand_t_cop[,3]) will give a confidence interval (95 percent by default) for the correlation (Pearson by default) between components 1 and 3. Does this confidence interval include 0.75?

   The question is asking why cor(ran_t_cop) is different from the correlation params passed into the copula cop_t_dim3.

   cor.test(rand_t_cop[,1],rand_t_cop[,3])

   ## 
   ##  Pearson's product-moment correlation
   ## 
   ## data:  rand_t_cop[, 1] and rand_t_cop[, 3]
   ## t = 22.314, df = 498, p-value < 2.2e-16
   ## alternative hypothesis: true correlation is not equal to 0
   ## 95 percent confidence interval:
   ##  0.6603249 0.7483624
   ## sample estimates:
   ##      cor 
   ## 0.707073

   The confidence interval is (0.6603, 0.7484) which does not include 0.75, which was passed into the copula object.

   0.75 is the correlation between the t-distributed random variables that define the copula and need not be the same as the uniformly-distributed variables in the copula itself.

Q3

# create Gaussian copula obj
cop_normal_dim3 = normalCopula(dim = 3, param = c(-0.6,0.75,0),
  dispstr = "un")
# define dist using mvdc (multivar dist constructed from copulas)
mvdc_normal = mvdc(copula = cop_normal_dim3, margins = rep("exp", 3),
  paramMargins = list(list(rate=2), list(rate=3), list(rate=4)))
# generate random sample
set.seed(5640)
rand_mvdc = rMvdc(n = 1000, mvdc = mvdc_normal)
pairs(rand_mvdc)

par(mfrow = c (2,2))
for(i in 1:3) plot(density(rand_mvdc[,i]))

1. What are the marginal distributions of the three components in rand mvdc? What are their expected values?

   The marginal distributions for all 3 components are Exponential with rate \(\lambda= 3\), which means that their mean value is \(1/\lambda = 1/3\).

2. Are the second and third components independent? Why or why not?

   The 2nd and 3rd components are (marginally) independent, since they jointly (i.e., var2 and var3 combined, without var1) follow 2D Normal/Gaussian, and their correlation is 0.

   Their correlation is given by the copula::p2P function, which specifies the correlation matrix based on the vector of distinct correlations; for this problem we had:

   # p2P() creates a matrix from a given vector of parameters. 
   # P2p() creates a numeric vector from a given matrix
   p2P( c( -0.6, +0.75, 0 ) )

   ##       [,1] [,2] [,3]
   ## [1,]  1.00 -0.6 0.75
   ## [2,] -0.60  1.0 0.00
   ## [3,]  0.75  0.0 1.00

   which means that

   • var1 & var2 are negatively correlated,
   • var1 & var3 are positively correlated, and
   • var2 & var3 are un-correlated

   Note that these relations are reflected in the Exponential random variates that were generated. Also note that the conditional distribution var2 and var3 given var1 would NOT be independent (that’s why we have to specify that the “marginal” distribution of var2 & var3 is independent.)

Q1

Write an R program to find the efficient frontier, the tangency portfolio, and the minimum variance portfolio, and plot on “reward-risk space” the location of each of the six stocks, the efficient frontier, the tangency portfolio, and the line of efficient portfolios.

Use the constraints that \(−0.1 \leq w_j \leq 0.5\) for each stock. The first constraint limits short sales but does not rule them out completely. The second constraint prohibits more than 50% of the investment in any single stock.

Assume that the annual risk-free rate is 3% and convert this to a daily rate by dividing by 365, since interest is earned on trading as well as nontrading days.

M = length(mean_vect)
library(quadprog)
Amat = cbind(rep(1,M),mean_vect,diag(1,nrow=M),-diag(1,nrow=M))
muP = seq(.05,0.08,length=300)
sdP = muP
weights = matrix(0,nrow=300,ncol=M)
# quadratic programming problem
for (i in 1:length(muP)){
  result = solve.QP(Dmat=cov_mat, dvec=rep(0,M), Amat=Amat,
    c(1,muP[i],rep(-.1,M),rep(-.5,M)), meq=2)
  sdP[i] = sqrt(2*result$value)
  weights[i,] = result$solution
}
plot(sdP,muP,type="l",xlim=c(0,2.5),ylim=c(0,.1))
mufree = 3/365
points(0,mufree,cex=3,col="blue",pch="*")
sharpe = (muP-mufree)/sdP
ind = (sharpe == max(sharpe)) # locates the tangency portfolio
weights[ind,] # weights of the tangency portfolio

## [1] -0.091181044 -0.002910879  0.335318542  0.383714329  0.319484849
## [6]  0.055574204

lines(c(0,sdP[ind]),c(mufree,muP[ind]),col="red",lwd=3)
points(sdP[ind],muP[ind],col="blue",cex=3,pch="*")
ind2 = (sdP == min(sdP))
points(sdP[ind2],muP[ind2],col="green",cex=3,pch="*")
ind3 = (muP > muP[ind2])
lines(sdP[ind3],muP[ind3],type="l",xlim=c(0,.25),
  ylim=c(0,.3),col="cyan",lwd=3)
text(sd_vect[1],mean_vect[1],"GM")
text(sd_vect[2],mean_vect[2],"F")
text(sd_vect[3],mean_vect[3],"UTX")
text(sd_vect[4],mean_vect[4],"CAT")
text(sd_vect[5],mean_vect[5],"MRK")
text(sd_vect[6],mean_vect[6],"IBM")
legend("topleft",c("efficient frontier","efficient portfolios",
  "tangency portfolio","min var portfolio"),
  lty=c(1,1,NA,NA),
  lwd=c(3,3,1,1),
  pch=c("","","*","*"),
  col=c("cyan","red","blue","green"),
  pt.cex=c(1,1,3,3)
  )

Q2

If an investor wants an efficient portfolio with an expected daily return of 0.07%, how should the investor allocate his or her capital to the six stocks and to the risk-free asset? Assume that the investor wishes to use the tangency portfolio computed with the constraints \(−0.1 \leq w_j \leq 0\)., not the unconstrained tangency portfolio.

Let \(\omega\) be the weight for the risk-free asset, \(\mu_f\) be the risk-free rate, and \(\mu_T\) be the expected return of the tangency portfolio. Then \(\omega\) solves \(0.07=\omega \mu_f+(1-\omega) \mu_T\). The following continuation of the \(\mathrm{R}\) program computes the weights for the six stocks and the risk-free asset. The last line checks that the weights sum to 1.

options(digits=3)
omega = (.07 - muP[ind]) / (mufree - muP[ind])
omega

## [1] 0.0518

1-omega

## [1] 0.948

(1-omega)*weights[ind]

## [1] -0.08645 -0.00276  0.31794  0.36382  0.30292  0.05269

omega + sum((1-omega)*weights[ind])

## [1] 1

Fitting Factor Models by Time Series Regression

In this section, we will start with the one-factor CAPM model of Chap. 17 and then extend this model to the three-factor Fama–French model. We will use the data set Stock_Bond_2004_to_2005.csv on the book’s website, which contains stock prices and other financial time series for the years 2004 and 2005. Data on the Fama–French factors are available at Prof. Kenneth French’s website where RF is the risk-free rate and Mkt.RF, SMB, and HML are the Fama–French factors. Go to Prof. French’s website and get the daily values of RF, Mkt.RF, SMB, and HML for the years 2004–2005. It is assumed here that you’ve put the data in a text file FamaFrenchDaily.txt. Returns on this website are expressed as percentages. Now fit the CAPM to the four stocks using the lm command. This code fits a linear regression model separately to the four responses. In each case, the independent variable is Mkt.RF.

start = as.Date("2004-01-01")
end = as.Date("2005-12-31")

stocks = read.csv("./res/datasets/Stock_Bond.csv",header=T)
stocks$Date <- as.Date(stocks$Date, format = "%d-%b-%y")
stocks_subset <- subset(stocks, Date >= start & Date <= end)
stocks_subset <- subset(stocks_subset, select = c("GM_AC", "F_AC", "UTX_AC", "MRK_AC"))

FF_data = read.table("./res/datasets/FamaFrenchDaily.txt", header = TRUE)
FF_data$date <- as.Date(as.character(FF_data$date), format = "%Y%m%d")
FF_data <- subset(FF_data, date >= start & date <= end)
FF_data = FF_data[-1, ] # delete first row since stocks_diff lost a row from differencing

stocks_diff = as.data.frame(100 * apply(log(stocks_subset),2, diff) - FF_data$RF)
names(stocks_diff) = c("GM", "Ford", "UTX", "Merck")
fit1 = lm(as.matrix(stocks_diff) ~ FF_data$Mkt.RF)
summary(fit1)

## Response GM :
## 
## Call:
## lm(formula = GM ~ FF_data$Mkt.RF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.832  -0.721   0.024   0.808  15.248 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -0.2290     0.0865   -2.65   0.0083 ** 
## FF_data$Mkt.RF   1.2500     0.1273    9.82   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.94 on 501 degrees of freedom
## Multiple R-squared:  0.161,  Adjusted R-squared:  0.16 
## F-statistic: 96.4 on 1 and 501 DF,  p-value: <2e-16
## 
## 
## Response Ford :
## 
## Call:
## lm(formula = Ford ~ FF_data$Mkt.RF)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.351 -0.851  0.018  0.790  9.182 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -0.1835     0.0676   -2.72   0.0069 ** 
## FF_data$Mkt.RF   1.3195     0.0995   13.26   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.51 on 501 degrees of freedom
## Multiple R-squared:  0.26,   Adjusted R-squared:  0.258 
## F-statistic:  176 on 1 and 501 DF,  p-value: <2e-16
## 
## 
## Response UTX :
## 
## Call:
## lm(formula = UTX ~ FF_data$Mkt.RF)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.353 -0.541 -0.007  0.537  3.613 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     0.00219    0.03883    0.06     0.95    
## FF_data$Mkt.RF  0.91932    0.05718   16.08   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.87 on 501 degrees of freedom
## Multiple R-squared:  0.34,   Adjusted R-squared:  0.339 
## F-statistic:  258 on 1 and 501 DF,  p-value: <2e-16
## 
## 
## Response Merck :
## 
## Call:
## lm(formula = Merck ~ FF_data$Mkt.RF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -31.155  -0.481   0.059   0.709  12.307 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -0.0888     0.0920   -0.97     0.33    
## FF_data$Mkt.RF   0.6255     0.1355    4.62    5e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.06 on 501 degrees of freedom
## Multiple R-squared:  0.0408, Adjusted R-squared:  0.0389 
## F-statistic: 21.3 on 1 and 501 DF,  p-value: 4.99e-06

Problem 5

The CAPM predicts that all four intercepts will be zero. For each stock, using α = 0.025, can you accept the null hypothesis that its intercept is zero? Why or why not? Include the p-values with your work.

The intercept’s p-value is below 0.025 for General motors (GM) and Ford, but not for United Technologies Incorporatted (UTX) and Merck.

Of course, a p-value only shows statistical significance, not the size of an effect. However, estimated intercepts for GM and Ford are −0.23 and −0.18 and these are reasonably large in magnitude. Since they are negative, this suggests than these two stocks were overpriced; see Section 17.6.3 for discussion.

Problem 6

The CAPM also predicts that the four sets of residuals will be uncorrelated. What is the correlation matrix of the residuals? Give a 95% confidence interval for each of the six correlations. Can you accept the hypothesis that all six correlations are zero?

cor(residuals(fit1))

##            GM     Ford     UTX    Merck
## GM     1.0000  0.52016 -0.0100 -0.08776
## Ford   0.5202  1.00000 -0.0238 -0.00958
## UTX   -0.0100 -0.02376  1.0000 -0.00550
## Merck -0.0878 -0.00958 -0.0055  1.00000

res = residuals(fit1)
cor.test(res[,"GM"],res[,"Ford"])

## 
##  Pearson's product-moment correlation
## 
## data:  res[, "GM"] and res[, "Ford"]
## t = 14, df = 501, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.453 0.581
## sample estimates:
##  cor 
## 0.52

cor.test(res[,"GM"],res[,"Merck"])

## 
##  Pearson's product-moment correlation
## 
## data:  res[, "GM"] and res[, "Merck"]
## t = -2, df = 501, p-value = 0.05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.17385 -0.00033
## sample estimates:
##     cor 
## -0.0878

All correlations are reasonably close to 0 (less than 0.1 in magnitude) except the correlation between GM and Ford residuals. That correlation is 0.52 and has a very small p-value.

The correlation between GM and Merck residuals is −0.0878 and is statistically significant at 0.05 but might be too small to be of practical significance.

Problem 7

Regardless of your answer to Problem 6, assume for now that the residuals are uncorrelated. Then use the CAPM to estimate the covariance matrix of the excess returns on the four stocks. Compare this estimate with the sample covariance matrix of the excess returns. Do you see any large discrepancies between the two estimates of the covariance matrix?

attach(FF_data)
sigF = var(Mkt.RF)
bbeta = as.matrix(fit1$coef)
bbeta = bbeta[-1,] # delete intercepts so bbeta has the four slopes
n=dim(stocks_diff)[1]
sigeps = as.matrix((var(as.matrix(res))))
sigeps = diag(as.matrix(sigeps))
sigeps = diag(sigeps,nrow=4)
cov_equities = sigF* bbeta %*% t(bbeta) + sigeps
cov_equities

##         GM  Ford   UTX Merck
## [1,] 4.464 0.761 0.530 0.360
## [2,] 0.761 3.090 0.559 0.381
## [3,] 0.530 0.559 1.145 0.265
## [4,] 0.360 0.381 0.265 4.423

cov(stocks_diff)

##           GM  Ford   UTX  Merck
## GM    4.4641 2.282 0.513 0.0108
## Ford  2.2825 3.090 0.528 0.3507
## UTX   0.5130 0.528 1.145 0.2553
## Merck 0.0108 0.351 0.255 4.4228

We see that the estimated covariance matrix using the CAPM is similar to the sample covariance matrix, with the exception of the covariance between GM and Ford. Since these two stocks have a high residual correlation and the CAPM assumes that the residual correlation is 0, it is not surprising that the CAPM estimated covariance matrix severely underestimates the correlation between GM and Ford.

Problem 8

Next, you will fit the Fama–French three-factor model. Run the following R code, which is much like the previous code except that the regression model has two additional predictor variables, SMB and HML.

fit2 = lm(as.matrix(stocks_diff) ~ FF_data$Mkt.RF + FF_data$SMB + FF_data$HML)
summary(fit2)

## Response GM :
## 
## Call:
## lm(formula = GM ~ FF_data$Mkt.RF + FF_data$SMB + FF_data$HML)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.759  -0.727  -0.033   0.775  15.007 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -0.2511     0.0865   -2.90   0.0038 ** 
## FF_data$Mkt.RF   1.3889     0.1517    9.16   <2e-16 ***
## FF_data$SMB     -0.2504     0.2236   -1.12   0.2633    
## FF_data$HML      0.6006     0.2643    2.27   0.0235 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.93 on 499 degrees of freedom
## Multiple R-squared:  0.173,  Adjusted R-squared:  0.168 
## F-statistic: 34.9 on 3 and 499 DF,  p-value: <2e-16
## 
## 
## Response Ford :
## 
## Call:
## lm(formula = Ford ~ FF_data$Mkt.RF + FF_data$SMB + FF_data$HML)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.344 -0.832 -0.010  0.793  9.198 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -0.1951     0.0679   -2.87   0.0042 ** 
## FF_data$Mkt.RF   1.3511     0.1191   11.35   <2e-16 ***
## FF_data$SMB     -0.0157     0.1756   -0.09   0.9288    
## FF_data$HML      0.3412     0.2074    1.64   0.1006    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.51 on 499 degrees of freedom
## Multiple R-squared:  0.264,  Adjusted R-squared:  0.26 
## F-statistic: 59.7 on 3 and 499 DF,  p-value: <2e-16
## 
## 
## Response UTX :
## 
## Call:
## lm(formula = UTX ~ FF_data$Mkt.RF + FF_data$SMB + FF_data$HML)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.269 -0.537  0.008  0.493  3.429 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -3.99e-06   3.88e-02    0.00   0.9999    
## FF_data$Mkt.RF  1.03e+00   6.80e-02   15.12   <2e-16 ***
## FF_data$SMB    -2.93e-01   1.00e-01   -2.92   0.0037 ** 
## FF_data$HML    -9.59e-04   1.19e-01   -0.01   0.9935    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.864 on 499 degrees of freedom
## Multiple R-squared:  0.352,  Adjusted R-squared:  0.348 
## F-statistic: 90.2 on 3 and 499 DF,  p-value: <2e-16
## 
## 
## Response Merck :
## 
## Call:
## lm(formula = Merck ~ FF_data$Mkt.RF + FF_data$SMB + FF_data$HML)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.627  -0.454   0.106   0.691  12.310 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -0.0598     0.0915   -0.65  0.51337    
## FF_data$Mkt.RF   0.7093     0.1605    4.42  1.2e-05 ***
## FF_data$SMB     -0.4174     0.2366   -1.76  0.07832 .  
## FF_data$HML     -0.9559     0.2796   -3.42  0.00068 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.04 on 499 degrees of freedom
## Multiple R-squared:  0.066,  Adjusted R-squared:  0.0604 
## F-statistic: 11.8 on 3 and 499 DF,  p-value: 1.86e-07

The CAPM predicts that for each stock, the slope (beta) for SMB and HML will be zero. Explain why the CAPM makes this prediction. Do you accept this null hypothesis? Why or why not?

Problem 9

If the Fama–French model explains all covariances between the returns, then the correlation matrix of the residuals should be diagonal. What is the estimated correlations matrix? Would you accept the hypothesis that the correlations are all zero?

Problem 10

Which model, CAPM or Fama–French, has the smaller value of AIC? Which has the smaller value of BIC? What do you conclude from this?

Problem 11

What is the covariance matrix of the three Fama–French factors?

Problem 12

In this problem, Stocks 1 and 2 are two stocks, not necessarily in the Stock_FX_Bond_2004_to_2005.csv data set. Suppose that Stock 1 has betas of 0.5, 0.4, and −0.1 with respect to the three factors in the Fama–French model and a residual variance of 23.0. Suppose also that Stock 2 has betas of 0.6, 0.15, and 0.7 with respect to the three factors and a residual variance of 37.0. Regardless of your answer to Problem 9, when doing this problem, assume that the three factors do account for all covariances.

1. Use the Fama–French model to estimate the variance of the excess return on Stock 1.

2. Use the Fama–French model to estimate the variance of the excess return on Stock 2.

3. Use the Fama–French model to estimate the covariance between the excess returns on Stock 1 and Stock 2.

Statistical Factor Models

This section applies statistical factor analysis to the log returns of 10 stocks in the data set Stock_FX_Bond.csv. The data set contains adjusted closing (AC) prices of the stocks, as well as daily volumes and other information that we will not use here. The following R code will read the data, compute the log returns, and fit a two-factor model. Note that factanal works with the correlation matrix or, equivalently, with standardized variables.

dat = read.csv("./res/datasets/Stock_FX_Bond.csv")
stocks_ac = dat[ , c(3, 5, 7, 9, 11, 13, 15, 17)]
n = length(stocks_ac[ , 1])
stocks_returns = log(stocks_ac[-1, ] / stocks_ac[-n, ])
fact = factanal(stocks_returns, factors = 2, rotation = "none")
print(fact)

## 
## Call:
## factanal(x = stocks_returns, factors = 2, rotation = "none")
## 
## Uniquenesses:
##   GM_AC    F_AC  UTX_AC  CAT_AC  MRK_AC  PFE_AC  IBM_AC MSFT_AC 
##   0.399   0.423   0.718   0.714   0.519   0.410   0.760   0.749 
## 
## Loadings:
##         Factor1 Factor2
## GM_AC    0.693  -0.348 
## F_AC     0.692  -0.313 
## UTX_AC   0.531         
## CAT_AC   0.529         
## MRK_AC   0.551   0.421 
## PFE_AC   0.574   0.511 
## IBM_AC   0.490         
## MSFT_AC  0.499         
## 
##                Factor1 Factor2
## SS loadings       2.64   0.666
## Proportion Var    0.33   0.083
## Cumulative Var    0.33   0.414
## 
## Test of the hypothesis that 2 factors are sufficient.
## The chi square statistic is 565 on 13 degrees of freedom.
## The p-value is 2.6e-112

Loadings less than the parameter cutoff are not printed. The default value of cutoff is 0.1, but you can change it as in “print(fact,cutoff = 0.01)” or “print(fact, cutoff = 0)”.

Problem 13

What are the factor loadings? What are the variances of the unique risks for Ford and General Motors?

We see from the output that the uniquenesses for Ford and GM are 0.423 and 0.399, respectively. The loadings are printed below (also shown in the output).

unique = as.numeric(fact$unique)
loadings = matrix(as.numeric(loadings(fact)), ncol = 2)
loadings

##       [,1]     [,2]
## [1,] 0.693 -0.34769
## [2,] 0.692 -0.31339
## [3,] 0.531 -0.02819
## [4,] 0.529 -0.07850
## [5,] 0.551  0.42090
## [6,] 0.574  0.51112
## [7,] 0.490  0.00424
## [8,] 0.499  0.04400

Problem 14

Does the likelihood ratio test suggest that two factors are enough? If not, what is the minimum number of factors that seems sufficient? The following code will extract the loadings and uniquenesses.

# try increasing the number of factors
factanal(stocks_returns, factors = 3, rotation = "none")

## 
## Call:
## factanal(x = stocks_returns, factors = 3, rotation = "none")
## 
## Uniquenesses:
##   GM_AC    F_AC  UTX_AC  CAT_AC  MRK_AC  PFE_AC  IBM_AC MSFT_AC 
##   0.402   0.368   0.713   0.703   0.510   0.386   0.577   0.622 
## 
## Loadings:
##         Factor1 Factor2 Factor3
## GM_AC    0.688  -0.339         
## F_AC     0.704  -0.334  -0.160 
## UTX_AC   0.528                 
## CAT_AC   0.528           0.112 
## MRK_AC   0.546   0.420  -0.126 
## PFE_AC   0.572   0.524  -0.111 
## IBM_AC   0.518           0.392 
## MSFT_AC  0.520           0.324 
## 
##                Factor1 Factor2 Factor3
## SS loadings      2.689   0.687   0.343
## Proportion Var   0.336   0.086   0.043
## Cumulative Var   0.336   0.422   0.465
## 
## Test of the hypothesis that 3 factors are sufficient.
## The chi square statistic is 162 on 7 degrees of freedom.
## The p-value is 1.06e-31

factanal(stocks_returns, factors = 4, rotation = "none")

## 
## Call:
## factanal(x = stocks_returns, factors = 4, rotation = "none")
## 
## Uniquenesses:
##   GM_AC    F_AC  UTX_AC  CAT_AC  MRK_AC  PFE_AC  IBM_AC MSFT_AC 
##   0.411   0.349   0.077   0.702   0.496   0.402   0.540   0.620 
## 
## Loadings:
##         Factor1 Factor2 Factor3 Factor4
## GM_AC    0.489   0.484  -0.333         
## F_AC     0.496   0.516  -0.345  -0.141 
## UTX_AC   0.936  -0.217                 
## CAT_AC   0.486   0.229                 
## MRK_AC   0.373   0.406   0.428  -0.129 
## PFE_AC   0.392   0.413   0.513  -0.107 
## IBM_AC   0.390   0.343           0.435 
## MSFT_AC  0.379   0.351           0.329 
## 
##                Factor1 Factor2 Factor3 Factor4
## SS loadings      2.185   1.176   0.685   0.356
## Proportion Var   0.273   0.147   0.086   0.044
## Cumulative Var   0.273   0.420   0.506   0.550
## 
## Test of the hypothesis that 4 factors are sufficient.
## The chi square statistic is 0.3 on 2 degrees of freedom.
## The p-value is 0.86

The results of the likelihood ratio tests below strongly suggest that there are more than two factors. It seems that four factors are sufficient, but not 3.

Problem 15

Regardless of your answer to Problem 14, use the two-factor model to estimate the correlation of the log returns for Ford and IBM

loadings = matrix(as.numeric(loadings(fact)),ncol=2)
unique = as.numeric(fact$unique)
options(digits=2)
loadings %*% t(loadings) + diag(unique)

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1.00 0.59 0.38 0.39 0.24 0.22 0.34 0.33
## [2,] 0.59 1.00 0.38 0.39 0.25 0.24 0.34 0.33
## [3,] 0.38 0.38 1.00 0.28 0.28 0.29 0.26 0.26
## [4,] 0.39 0.39 0.28 1.00 0.26 0.26 0.26 0.26
## [5,] 0.24 0.25 0.28 0.26 1.00 0.53 0.27 0.29
## [6,] 0.22 0.24 0.29 0.26 0.53 1.00 0.28 0.31
## [7,] 0.34 0.34 0.26 0.26 0.27 0.28 1.00 0.24
## [8,] 0.33 0.33 0.26 0.26 0.29 0.31 0.24 1.00

cor(stocks_returns)

##         GM_AC F_AC UTX_AC CAT_AC MRK_AC PFE_AC IBM_AC MSFT_AC
## GM_AC    1.00 0.62   0.35   0.36   0.25   0.23   0.32    0.31
## F_AC     0.62 1.00   0.35   0.37   0.26   0.25   0.30    0.30
## UTX_AC   0.35 0.35   1.00   0.40   0.26   0.28   0.29    0.28
## CAT_AC   0.36 0.37   0.40   1.00   0.24   0.25   0.30    0.29
## MRK_AC   0.25 0.26   0.26   0.24   1.00   0.55   0.24    0.27
## PFE_AC   0.23 0.25   0.28   0.25   0.55   1.00   0.26    0.29
## IBM_AC   0.32 0.30   0.29   0.30   0.24   0.26   1.00    0.41
## MSFT_AC  0.31 0.30   0.28   0.29   0.27   0.29   0.41    1.00

(a)

Download the adjusted daily closing prices from Jan 1 2020 to Dec 31 2022 for the 5 mutual funds below (use tseries::get.hist.quote() for each ticker):

• FCNTX: Fidelity Contrafund
• PIODX: Pioneer A
• AIVSX: American Funds Invmt Co of Amer A
• PRBLX: Parnassus Core Equity Investor
• VFIAX: Vanguard 500 Index Admiral

Note that each of these funds has at least 90% of their weight in the US stocks market. You can actually check the composition of the investment over different stock sectors from Yahoo Finance, under the fund’s holdings tab; e.g. for FCNTX at https://finance.yahoo.com/quote/FCNTX/holdings.

library(zoo)
library(tseries)

MF.names=c('FCNTX','PIODX','AIVSX','PRBLX','VFIAX')
N.MF=length(MF.names)
S=list()
for(i in 1:N.MF){
  S[[i]]=get.hist.quote(MF.names[i], start='2022-01-01',
                        end='2022-12-31', quote='AdjClose', quiet = TRUE)
}
R=lapply(S, FUN = function(x){ diff(x) / timeSeries::lag(x,-1) }) # calculate MF returns
RY=matrix(unlist(R),ncol=N.MF) # bind returns in a matrix
colnames(RY)=MF.names

(b)

Assume you do not have any information about the investment strategy of the funds.Download the daily prices and calculate returns of the following EFTs, which track different sectors of the economy:
• XLB: Basic Materials
• XLY: Consumer Cyclical
• XLF: Financial Services
• VNQ: Real Estate
• XLP: Consumer Defensive
• XLV: Healthcare
• XLU: Utilities
• XTL: Communication Services
• XLE: Energy
• XLI: Industrials
• XLK: Technology

Regress each of the mutual fund returns on the above ETF returns and create barplots of the estimated beta coefficients. Do these accurately reflect the allocation over the different sectors (as described in Yahoo Finance)?

ETF.names=c('XLB','XLY','XLF','VNQ','XLP','XLV','XLU','XTL','XLE','XLI','XLK')
sectors = c("Basic Materials","Consumer Cyclical","Financial Services",
            "Real Estate","Consumer Defensive", "Healthcare", "Utilities", 
            "Communication Services", "Energy", "Industrials", "Technology")
N.ETF=length(ETF.names)
S=list()
for(i in 1:N.ETF){
  S[[i]]=get.hist.quote(ETF.names[i], start='2022-01-01',
                        end='2022-12-31', quote='AdjClose', quiet = TRUE)
}

## Warning: XLE download failed; trying again.

R=lapply(S, FUN = function(x){ diff(x) / timeSeries::lag(x,-1) }) # calculate ETF returns
RX=matrix(unlist(R),ncol=N.ETF) # bind returns in a matrix
colnames(RX)=ETF.names

out=list(); 
for(i in 1:N.MF){
  out[[i]]=lm(RY[,i] ~ RX) # regress mutual funds returns on ETF returns
  weights = out[[i]]$coef[-1] 
  barplot( rev(weights), names.arg = rev(sectors), main = MF.names[i],
          horiz = TRUE, las = 2, cex.names = .8)
}

The barplots of the regression coefficients (betas) roughly follow the sector weightings for each fund. Nevertheless, they are not always close in actual value (e.g. in some cases the betas are negative, even though weightings are positive). The differences can be due to the fact that we use ETFs as proxies for a sector, but the actual holding of the fund within the sector might be different. Moreover, there will be estimation error in our regression model, which is only based on the last year’s returns.

Nevertheless, it is quite impressive that we can (approximately) identify the strategy of a fund, without knowing anything beyond its past returns. This approach works because of the linear formula for net portfolio returns: \[ R_p = w_1 R_1 + \cdots + w_N R_N\] Regressing portfolio returns on other assets, we can estimate the weights, assuming the portfolio composition is constant.

(c)

Compare the performance of the mutual funds to that of a portfolio of ETFs by reporting the value of Jensen’s alpha (based on the regressions from the previous part) and its corresponding p-value.

alpha=p.val=rep(0,N.MF)
for(i in 1:N.MF){
  alpha[i]=out[[i]]$coefficients[1]*250 # annualized Jensen alpha
  p.val[i]=summary(out[[i]])$coefficients[1,4]
}
cbind(alpha, p.val)

##       alpha p.val
## [1,] -0.103 0.078
## [2,] -0.034 0.370
## [3,] -0.039 0.250
## [4,] -0.025 0.468
## [5,] -0.016 0.348

All the funds’ alphas are negative, although their p-values are not very small. A likely cause for this is that funds charge a fee which consistently eats up some of the returns of their constituent assets. ETFs typically have lower fees than mutual funds, but our regression does not account for transaction costs (it is more costly to buy multiple assets than a single one), so the comparison is more nuanced.

Note that you can find Jensen alphas and other performance measures (e.g., Sharpe & Treynor Ratios) for assets in Yahoo! Finance under the risk tab. These metrics are based on the CAPM/Market factor model, by regressing the asset’s returns on a proxy for the market return (e.g., S&P500).

(a) [10 points]

Simulate \(n=10,000\) random variates from a \(\operatorname{GPD}(\gamma=.5, \sigma=1)\) distribution using the inverse CDF method from Q1.(d), and plot the histogram of their log-values.?

set.seed(123)
n = 10000
gamma = .5
sigma = 1

X = sigma/gamma * ( (1-runif(n))^(-gamma) - 1 ) 
hist( log( X ) )

(b) [15 points]

Simulate another 10, 000 random variates from \(\operatorname{GPD}(\gamma=.5, \sigma=1)\), but now use the mixture method from Q1.(e). Create a QQ-plot of the two sets of variates (from this and the previous part0), and comment whether they seem to come from the same distribution? (Hint: Use qqplot \((\mathrm{X}, \mathrm{Y})\), where \(\mathrm{X}, \mathrm{Y}\) are the values you generated.)

alpha = 1 / gamma
beta = sigma / gamma

L = rgamma( n, shape = alpha, rate = beta )
Y = rexp( n, rate = L)

qqplot(X,Y); abline( c(0,1) )

(c) [15 points]

Take the values from Q2.(a) and calculate their exceedances above 10, i.e. \(Z=X-10\) for values of \(X\) that are greater than 10. From Q1.(b), we know that \(Z\) must follow \(G P D(\gamma, \sigma+10 \gamma)\). Verify that the exceedances follow this distribution, by creating a QQ-plot of \(Z\)-values versus their theoretical quantiles.

u = 10
# theoretical quantiles
Z = (X-u)[X>u] # values of X-u that satisfy X>u
Zs = sort(Z) 
# sample quantiles
nZ = length(Z)
sigma_u = sigma + gamma * u
Q = sigma_u/gamma * ( (1 - ppoints(nZ) )^(-gamma) - 1 ) 
# ppoints(n) creates n increments of size 1/n
plot( Zs, Q ); abline(0,1)

(d) [10 points]

Finally, verify through simulation that conditional exceedances of fat tail distributions with tail index \(\alpha\) approach a GPD with \(\gamma=1 / \alpha\). Generate \(n=10,000\) values from a \(t(d f=4)\) distribution (i.e. \(\alpha=4\) ) and take their absolute value (to avoid wasting values by symmetry). Calculate the exceedances above 10 again, and create their QQ-plot versus quantiles from a \(\operatorname{GPD}(\gamma=1 / \alpha, \sigma=1)\). (Note: the points in your QQ-plot should lie close to a straight line, but the slope does not need to be equal to 1 because of the arbitrary GPD scale \(\sigma=1\) ).

# theoretical quantiles
W = abs( rt( 100000, 4 ) )
Ws = (W-u)[W>u]
Ws = sort(Ws)
# sample quantiles
nW = length(W)
sigma = 1
gamma = 1/4
Q = (1+u) * ( (1 - ppoints(nW))^(-gamma) - 1 ) # why (1+u)?
# Q = ((sigma + gamma * u) / gamma) * ( (1 - ppoints(nW))^(-gamma) - 1 )
qqplot( Ws, Q ); abline(0, 1)

part (a)

[4 points] Simulate 10, 000 Cauchy random variates and plot their cumulative mean versus the sample size, as in the plot below

n = 10000
set.seed(123)
x = rcauchy(n)
m = cumsum(x)/(1:n)
plot(m, type = 'l' )

part (b)

[8 points] Repeat the following experiment N = 1000 times: simulate n = 100 Cauchy variates and calculate their mean (in the end you should have 1,000).

From question 1, we know these means should be also Cauchy distributed. Verify this by creating a QQ-plot of the simulated values versus their theoretical quantiles, and comment on the fit.

Hint: Use qcauchy and ppoints() to generate theoretical quantiles from the Cauchy distribution.)

N = 1000; n = 100
x = replicate( n = N, expr = mean(rcauchy(n = n)) )
qqplot( sort(x), qcauchy( p = ppoints(N)) ); abline(c(0,1))

part (c)

[8 points] Now consider the Pareto distribution, with CDF \(F_X(x)=1-(1 / x)^\alpha, \forall x>1\) Use the inverse CDF method to generate 10,000 Pareto random variates with tail index \(\alpha=3\), and plot their cumulative mean vs sample size, as in part (a).

n = 10000
alpha = 3
u = runif(n)
y = (1 - u)^(-1/alpha)
my = cumsum(y)/(1:n)
plot(my, type = 'l' )
abline(h=1.5, col = 2, lty=2)

part (a)

[5 points] Calculate the variance-covariance matrix of the returns \(\boldsymbol{\Sigma}_R=\mathbb{V}[\boldsymbol{R}]\).

beta = matrix( c(1,2,3,3,2,1)/10, nrow = 2, byrow = TRUE )
S_e = diag( c(.1,.2) )
S_R = beta %*% t(beta) + S_e

part (b)

[5 points] Calculate the mean and variance of the return of the minimum-variance portfolio.

Si = solve(S_R) # Sigma inverse
w = rowSums(Si)/sum(Si)
mu = c(.03, .04)

mvp_mu = sum(w*mu)
mvp_s =  t(w) %*% S_R %*% w # also equal to 1/sum(Si)

part (c)

[10 points] Assume the risk-free rate is .015, and find the mean and variance of the tangency portfolio.

rf = .015
w_tp = Si %*% (mu - rf) / sum( Si %*% (mu - rf) )

sum( w_tp * mu ) # tangency portfolio mean 

## [1] 0.03634

t(w_tp) %*% S_R %*% w_tp # tangency portfolio variance 

##        [,1]
## [1,] 0.2152

part (d)

[10 points] Assume you could eliminate one of the factors (i.e. set \(\mathbb{V}\left[F_i\right]=0\) ). Which factor would you choose to eliminate in order to further reduce the variance of the minimumvariance portfolio? (justify your answer with a calculation & comparison)

mvp_vars = rep(0,3)

for(i in 1:3){
  beta_t = beta[,-i] # t for two factor
  S_t = beta_t %*% t(beta_t) + S_e
  Si_t = solve(S_t)
  mvp_vars[i] = 1/sum(Si_t) # min variance
}

which.min(mvp_vars)

## [1] 3

(a)

Simulate \(n=1,000\) paths of \(m=252\) daily returns each, from an AR(1) model with \(R_0=0, \varphi=.25, \mu=.05 / 252, \sigma_{\varepsilon}=.2 / \sqrt{252}, r_f=.02 / 252\).

n=1000; m=252
phi=.25; mu=.05/252; sigma=.4/sqrt(252); rf=.02/252
V0=1

set.seed(1234567890)
R=matrix(0,m,n)
Z=matrix(rnorm(n*m),m,n)
# Initialize returns at stationary distribution
R[1,]=mu+(0-mu)*phi+sigma*Z[1,] 
for(i in 2:m){
  R[i,]=mu+(R[i-1,]-mu)*phi+sigma*Z[i,]
}
# To avoid the loop (& speed up computation) you can use the filter() function: 
# R = apply( mu*(1-phi) + sigma*Z, 2, filter, filter=phi, method="recursive")

(b)

Apply the static Kelly criterion based on the stationary mean \(\mathbb{E}\left[R_t\right]=\mu\) and variance \(\mathbb{V}\left[R_t\right]=\frac{\sigma_{\varepsilon}^2}{\left(1-\varphi^2\right)}\) of the \(\mathrm{AR}(1)\) process, and calculate the resulting wealth process for each path. Plot the wealth of the strategy on the first path, and the histogram of the final log-wealths across all paths.

sigma.b=sigma/sqrt(1-phi^2)
(f.b=(1+rf)*(mu-rf)/(sigma.b^2)) # static Kelly fraction

## [1] 0.1758

V.b=matrix(0,m,n)
V.b[1,]=V0*(1+rf+f.b*(R[1,]-rf))
for(i in 2:m){
  V.b[i,]=V.b[i-1,]*(1+rf+f.b*(R[i,]-rf))
}
# To avoid the loop use
# V.b=apply( 1+rf+f.b*(R-rf), 2, cumprod)

plot(V.b[,1], type='l')

hist( V.b[m,], main = "Static Kelly", xlab="Vn" )

Using the “static” Kelly criterion, we do not really take advantage of the predictability of the process, just of the fact that the asset returns are on average better than the risk-free rate. The average final wealth over 1 year is a modest 1.0247.

(c)

Apply the dynamic Kelly criterion and calculate the resulting wealth process for each path. Plot the first path and the histogram of the final log-wealths across all paths, and compare them to those of the previous part.

V.c=matrix(1,m,n)

# Dynamic Kelly fraction
mu_i=(0-mu)*phi 
(f.c_i=(1+rf)*mu_i/(sigma^2)) # if f<0 -> short asset 

## [1] -0.07813

V.c[1,]=V0*(1+rf+f.c_i*(R[1,]-rf))
for(i in 2:m){
  # Dynamic Kelly fraction
  mu_i=(R[i-1,]-mu)*phi 
  f.c_i=(1+rf)*mu_i/(sigma^2) # changing fraction 
  V.c[i,]=V.c[i-1,]*(1+rf+f.c_i*(R[i,]-rf))
}
# To avoid the loop use
# mu_i=( rbind(rep(0,n),R[-m,])-mu)*phi
# f.c_i=(1+rf)*mu_i/(sigma^2+mu_i^2) 
# V.c=apply( 1+rf+f.c_i*(R-rf), 2, cumprod)

# Plot wealth of 1st path 
plot(V.c[,1], type='l', main = "Wealth path")
lines(V.b[,1], type='l', col=2)
legend( "topleft", lwd = 2, col = 1:2, c("dynamic","static"))

# Plot log-wealth, for better comparison. 
plot(V.c[,1], type='l', log = "y", main = "log-Wealth path")
lines(V.b[,1], type='l', col=2)
legend( "topleft", lwd = 2, col = 1:2, c("dynamic","static"))

The dynamic Kelly wealth is growing much faster than the static one. That’s because it can take advantage of the predictability of the price on each day, and even when the price goes down, by using a “negative” fraction (i.e. shorting the asset). However, it is possible for the dynamic Kelly wealth to become negative if the optimal fraction is >1 (i.e., you use leverage). In fact, 22.4% of the dynamic Kelly paths had negative wealth (i.e. ruin/bankruptcy) at some point. But for those paths that did not go negative, the dynamic strategy led to much higher terminal wealth.

par(mfrow=c(1,2))
hist(log(V.b[m,]), main = "Static Kelly", xlab = "log(Vn)" ) 
hist(log(V.c[m, V.c[m,]>0]), main = "Dynamic Kelly", xlab = "log(Vn[Vn>0])" ) 

One can moderate and improve the dynamic strategy by constraining the invested fraction to be always \(|f^*|<1\) (i.e. avoid leverage). This technically removes the risk of ruin while still giving superior results:

V.c.constr=matrix(1,m,n)
mu_i=(0-mu)*phi 
f.c_i= min( max( (1+rf)*mu_i/(sigma^2), -1 ), 1) 
V.c[1,]=V0*(1+rf+f.c_i*(R[1,]-rf))
for(i in 2:m){
  # Dynamic Kelly fraction
  mu_i=(R[i-1,]-mu)*phi 
  f.c_i= pmax( pmin( (1+rf)*mu_i/(sigma^2), 1) , -1) 
  V.c.constr[i,]=V.c.constr[i-1,]*(1+rf+f.c_i*(R[i,]-rf))
}

par(mfrow=c(1,2))
hist(log(V.b[m,]), main = "Static Kelly", xlab = "log(Vn)" ) 
hist(log(V.c.constr[m,]), main = "Dynamic Kelly", xlab = "log(Vn)" ) 

(a) floating strike

Perform Monte Carlo simulation for estimating the price of a floating strike European lookback call option. Assume the underlying asset price follows Geometric BM with \(S_0=90, \quad T=1, r=2 \%, \sigma=20 \%\), and use unbiased estimation (i.e. simulate from the exact distribution of the minimum). Use \(n=100,000\) samples and create a \(95 \%\) confidence interval. Compare your result with the exact price of the option, which is 14.26674

set.seed(12345)

n=1000000 # number of paths
r=.02 # risk-free rate
v=.2 # volatility
M=1 # maturity
S0=90 # starting asset price 

# use reflection trick for simulating minimum:
# minimum of arithmetic BM is minus the maximum of process with opposite drift

Z=rnorm(n)
XT.reflected=-(r-v^2/2)*M + v*sqrt(M)*Z # use negative drift for reflected path
MT.reflected=(XT.reflected + sqrt( XT.reflected^2 -2*v^2*M*log(runif(n))))/2 # reflected max
ST=S0*exp(-XT.reflected) # final value
mT=S0*exp(-MT.reflected) # minimum 
payoff_2.a = exp(-r*M) * (ST-mT)
mean_2.a=mean(payoff_2.a) # MC estimate
se_2.a=sd(payoff_2.a)/sqrt(n) # st. dev
CI_2.a=mean_2.a+c(-1,1)*qnorm(.975)*se_2.a # 95% CI

print( paste("MC Float Strike Lookback call price =", mean_2.a) )

## [1] "MC Float Strike Lookback call price = 14.2733852982329"

print( paste( c("MC Float Strike Lookback call 95% CI =", CI_2.a), collapse = "  " ) )

## [1] "MC Float Strike Lookback call 95% CI =  14.2487087208825  14.2980618755833"

print( paste("Exact Float Strike Lookback call price =", 14.2667437) )

## [1] "Exact Float Strike Lookback call price = 14.2667437"

(b) fixed strike

Perform Monte Carlo simulation for estimating the price of a fixed strike European lookback call option. Assume the underlying asset price follows Geometric BM with \(S_0=90, \quad K=90, \quad T=1, r=2 \%, \sigma=20 \%\), and use unbiased estimation (i.e. simulate from the exact distribution of the maximum). Use \(n=100,000\) samples and create a \(95 \%\) confidence interval. Compare you result with the exact price of the option, which is 16.04886

K=90 # fixed strike
Z=rnorm(n)
XT=(r-v^2/2)*M + v*sqrt(M)*Z 
MT=S0 * exp ( (XT+sqrt(XT^2 -2*v^2*M*log(runif(n))))/2 )
payoff_2.b = exp(-r*M)*pmax(MT-K,0)
mean_2.b=mean(payoff_2.b) # MC estimate
se_2.b=sd(payoff_2.b)/sqrt(n) # st. dev
CI_2.b=mean_2.b+c(-1,1)*qnorm(.975)*se_2.b # 95% CI

print( paste("MC Fixed Strike Lookback call price =", mean_2.b) )

## [1] "MC Fixed Strike Lookback call price = 16.0314639209374"

print( paste( c( "MC Fixed Strike Lookback call 95% CI =", CI_2.b), collapse = "  " ) )

## [1] "MC Fixed Strike Lookback call 95% CI =  16.0052209559405  16.0577068859343"

print( paste("Exact Fixed Strike Lookback call price =", 16.0488631) )

## [1] "Exact Fixed Strike Lookback call price = 16.0488631"

(c) Euler discretization

Repeat the previous part assuming the risk-neutral asset price dynamics: \[ d S_t=\mu\left(t, S_t\right) d t+\sigma\left(t, S_t\right) d W_t=\left(r S_t\right) d t+\left(\sigma \cos \left(e^t\right) \sqrt{S_t}\right) d W_t \] Note that the price process is not a Geometric Brownian Motion anymore. Use Euler discretization with \(m=50\) steps and simulate \(n=10,000\) paths. In order to approximate the maximum of each path, simulate the maximum of each step using the result for the extrema of an arithmetic BM: \[ M_{t_i} \mid S_{t_{i-1}}, S_{t_i}=S_{t_{i-1}}+\frac{\Delta S_{t_i}+\sqrt{\Delta S_{t_i}^2-2 \cdot \Delta t \cdot \sigma^2\left(t_{i-1}, S_{t_{i-1}}\right) \cdot \log (U)}}{2} \] Where \[ M_{t_i}=\max _{t_{i-1} \leqslant t \leqslant t_i}\left\{S_t\right\}, \Delta S_{t_i}=S_{t_i}-S_{t_{i-1}}, \Delta t=t_i-t_{i-1} \& U \sim \text { Uniform }(0,1) \] The maximum of the entire path will be the maximum over all steps in the path, i.e. \(M_T=\max _{0 \leqslant t \leqslant T}\left\{S_t\right\}=\max _{i=0, \ldots, m}\left\{M_{t_i}\right\}\). Include a \(95 \%\) confidence interval with your answer.

K=90
n=10000; m=50; Dt=M/m
ST=matrix(S0,n,m+1); MT=matrix(0,n,m)
for(i in 1:m){
  ST.prev=ST[,i]
  step.mu=r*ST.prev*Dt; 
  step.sig=v*cos(exp(i*Dt))*sqrt(ST.prev)*sqrt(Dt)
  DST=step.mu+step.sig*rnorm(n)
  ST[,i+1]=ST[,i]+DST
  MT[,i]=ST[,i]+(DST+sqrt(DST^2-2*step.sig^2*log(runif(n))))/2
}
MT.all=apply(MT,1,max)
ST.final=ST[,m+1]
payoff_2.c = exp(-r*M)*pmax(MT.all-K,0)
mean_2.c=mean(payoff_2.c) # MC estimate
se_2.c=sd(payoff_2.c)/sqrt(n) # st. dev
CI_2.c=mean_2.c+c(-1,1)*qnorm(.975)*se_2.c # 95% CI

print( paste("MC Fixed Strike Lookback call price =", mean_2.c) )

## [1] "MC Fixed Strike Lookback call price = 2.18656942157468"

print( paste( c( "MC Fixed Strike Lookback call 95% CI =", CI_2.c), collapse = "  " ) )

## [1] "MC Fixed Strike Lookback call 95% CI =  2.17236758314114  2.20077126000821"

(a) GBM price

Estimate the price of a European gap call option with \(S_0=55, K_1=50, K_2=60, T=1 y r, r=2 \%, \sigma=20 \%\) and assuming the price follows GBM. Use \(n=100,000\) paths and calculate the standard error of your result. Compare that to the exact price of 6.129964.

n=100000; S0=55; K1=50; K2=60; M=1; r=.02; v=.2

Z=rnorm(n)
ST=S0*exp( (r-v^2/2)*M + v*sqrt(M)*Z)
payoff_2.a=exp(-r*M)*(ST-K1)*(ST>K2) 

(mean_2.a=mean(payoff_2.a))

## [1] 6.121

(se_2.a=sd(payoff_2.a)/sqrt(n))

## [1] 0.03081

(b) antithetic variables

Repeat part (a) using antithetic variables, and report the new price estimate and standard error.

Z=rnorm(n/2); Z=cbind(Z,-Z)
ST=S0*exp( (r-v^2/2)*M + v*sqrt(M)*Z)
payoff_2.b=rowMeans(exp(-r*M)*(ST-K1)*(ST>K2))
(mean_2.b=mean(payoff_2.b))

## [1] 6.121

(se_2.b=sd(payoff_2.b)/sqrt(n))

## [1] 0.01694

(c) stratification

Repeat part (a) using stratification with m = 5, 10, 20 equiprobable strata, and report the new price estimate and standard error

m.all=c(5,10,20)
mean_2.c=se_2.c=rep(0,length(m.all))
# iterate over 3 strata values
for(i in 1:length(m.all)){
  m=m.all[i]
  nj=n/m # num paths / num strata
  Z=matrix(0,nj,m)
  # iterate over each of m strata
  for(j in 1:m){
    temp=runif(nj,(j-1)/m,j/m)
    Z[,j]=qnorm(temp)
  } 
  ST=S0*exp( (r-v^2/2)*M + v*sqrt(M)*Z )
  payoff_2.c=exp(-r*M)*(ST-K1)*(ST>K2)
  mean_2.c[i]=mean( colMeans(payoff_2.c) )
  se_2.c[i]=sqrt( sum( apply(payoff_2.c,2,var)/nj ) )/m 
}
cbind(mean_2.c, se_2.c)

##      mean_2.c   se_2.c
## [1,]    6.124 0.013190
## [2,]    6.129 0.008517
## [3,]    6.131 0.005919

(d) comparison

Would you expect the variance reduction methods in part (b) and (c) to work better for in-the-money, at-the-money, or out-of-the-money options? Justify your answer

Antithetic variables: works best when original and antithetic variates are negatively related

Stratification: works best for highly variable payoffs (if target RV changes over its domain)

(a) GBM price

Estimate the price of an at-the-money forward start European call option with \(T_{\exp }=2 y r, r=2 \%, \sigma=20 \%\) and assuming the price follows GBM. Use \(n=10,000\) paths and calculate the standard error of your result. You should use the BlackScholes formula to explicitly calculate the value of the option at \(T_{\text {start }}=1\), i.e. you should generate \(S_1\) and then calculate its “payoff” as the Black-Scholes price of a European at-the-money call option with 1 year to maturity (use the BSprice() function provided in EX9.R). Compare that to the exact price of 6.241225.

n=10000; S0=70; Tstart=1; Texp=2; r=.02; v=.2

Z=rnorm(n)
STstart=S0*exp( (r-v^2/2)*Tstart + v*sqrt(Tstart)*Z )
payoff_3.a = exp(-r*Tstart) * BSprice(STstart, STstart, r, Texp-Tstart, v)[,1]
(mean_3.a=mean(payoff_3.a)) # MC estimate

## [1] 6.246

(se_3.a=sd(payoff_3.a)/sqrt(n))

## [1] 0.01266

(b) antithetic variables

Repeat part (a) using antithetic variables, and report the new price estimate and standard error.

Z=rnorm(n/2); Z.a=-Z
STstart=S0*exp( (r-v^2/2)*Tstart + v*sqrt(Tstart)*Z )
STstart.a=S0*exp( (r-v^2/2)*Tstart + v*sqrt(Tstart)*Z.a ) # antithetic path
payoff=exp(-r*Tstart) * BSprice(STstart, STstart, r, Texp-Tstart, v)[,1]
payoff.a=exp(-r*Tstart) * BSprice(STstart.a, STstart.a, r, Texp-Tstart, v)[,1] # antithetic path payoff
payoff_3.b=(payoff+payoff.a)/2
(mean_3.b=mean(payoff_3.b)) # MC estimate

## [1] 6.244

(se_3.b=sd(payoff_3.b)/sqrt(length(payoff_3.b))) # std. error

## [1] 0.002629

(c) stratification

Repeat part (a) using stratification with equiprobable strata, and report the new price estimate and standard error.

m=10; nj=n/m
Z=STstart=payoff_3.c=matrix(0,nj,m)
for(j in 1:m){
  temp=runif(nj,(j-1)/m,j/m)
  Z[,j]=qnorm(temp)
  STstart[,j]=S0*exp( (r-v^2/2)*Tstart + v*sqrt(Tstart)*Z[,j] )
  payoff_3.c[,j]=exp(-r*Tstart) * BSprice(STstart[,j], STstart[,j], r, Texp-Tstart, v)[,1]
} 
(mean_3.c=mean( colMeans(payoff_3.c) ) )

## [1] 6.239

(se_3.c=sqrt( sum( apply(payoff_3.c,2,var)/nj ) )/m )

## [1] 0.002762

(d) proof

(paper & pencil part) Show that the price of the at-the-money forward call option is actually a simple linear function of \(S_0\)

(Hint: use risk-neutral pricing and write the payoff of the option at \(T_{\text {start }}\) as the Black-Scholes price of an at-the-money call with maturity at \(T>T_{\text {start }}\) )

Dt=Texp-Tstart
d1=(r/v+v/2)*sqrt(Dt)
d2=(r/v-v/2)*sqrt(Dt)
(exact_price=S0*( pnorm(d1) - exp(-r*Dt)*pnorm(d2) ))

## [1] 6.241

The Black-Scholes price for a European call is \[ C=S_0 \Phi\left(d_1\right)-e^{-r T} K \Phi\left(d_2\right) \] where \(d_1=\frac{\ln \left(S_0 / K\right)+\left(r+\frac{1}{2} \sigma^2\right) T}{\sigma \sqrt{T}} ; d_2=d_1-\sigma \sqrt{T}\).

For an at-the-money call, in particular, we get \(C=S_0\left[\Phi\left(d_1\right)-e^{-r T} \Phi\left(d_2\right)\right]=\) \(S_0\left[\Phi\left(\frac{\left(r+\frac{1}{2} \sigma^2\right) T}{\sigma \sqrt{T}}\right)-e^{-r T} \Phi\left(\frac{\left(r-\frac{1}{2} \sigma^2\right) T}{\sigma \sqrt{T}}\right)\right]\)

By risk-neutral pricing, the price of the at-the-money forward call option becomes: \[ \begin{aligned} C_{F W D} & =\mathbb{E}\left[e^{-r T_{\exp }} \text { payoff }\right] \\ & =\mathbb{E}\left\{\mathbb{E}\left[e^{-r T_{\exp }}\left(\text { at the money Euro call at } T_{\text {start }}\right) \mid S_{T_{\text {start }}}\right]\right\} \\ & =\mathbb{E}\left[e^{-r T_{\text {start }}} C\left(K=S_{T_{\text {start }}}, T=T_{\exp }-T_{\text {start }}\right)\right]\\ & =\alpha \mathbb{E}\left[e^{-r T_{\text {start }}} S_{T_{\text {start }}}\right]=\alpha S_0 \end{aligned} \] where \[ \alpha = \Phi\left(\frac{\left(r+\frac{1}{2} \sigma^2\right)\left(T_{\exp }-T_{\text {start }}\right)}{\sigma \sqrt{T_{\exp }-T_{\text {start }}}}\right)-e^{-r\left(T_{\exp }-T_{\text {start }}\right)} \Phi\left(\frac{\left(r-\frac{1}{2} \sigma^2\right)\left(T_{\exp }-T_{\text {start }}\right)}{\sigma \sqrt{T_{\exp }-T_{\text {start }}}}\right) \] In the last line we used the fact that the discounted asset price is a martingale under the risk-neutral measure.

(a) path discretization

Perform simple Monte Carlo using path discretization with \(t_j=j \frac{T}{m}, j=0, \ldots, m\) and \(m=50\) in order to approximate the continuous median \(\left\{S_t\right\}_{t=0}^T\) by the discrete median \(\left\{S\left(t_j\right)\right\}_{j=0}^m\). Use \(n=10,000\) paths and report the standard error of the estimate.

S0=60; K=70; M=1; r=.05; v=.25
n=10000; m=50; Dt=M/m

Z=matrix(rnorm(n*m), n, m)
S=S0*exp(t(apply((r-v^2/2)*Dt+Z*v*sqrt(Dt),1,cumsum))) 
medS=apply(S,1,median)
payoff_4.a=exp(-r*M)*pmax(medS-K,0)
(mean_4.a=mean(payoff_4.a))

## [1] 1.036

(se_4.a=sd(payoff_4.a)/sqrt(n))

## [1] 0.03346

(b) control variates

Repeat part (a) using the Black-Scholes price of the European call as a control variate, and report the new price estimate and standard error

EX=BSprice(S0, K, r, M, v)[1]
X=exp(-r*M)*pmax(S[,m]-K,0) # control variate values
Y=payoff_4.a
(rho_XY=cor(X,Y))

## [1] 0.7177

(b=sd(Y)/sd(X)*rho_XY)

## [1] 0.3005

(mean_4.b=mean_4.a - b*(mean(X)-EX))

## [1] 1.039

(se_4.b=se_4.a*sqrt((1-rho_XY^2)))

## [1] 0.0233

(a) path discretization

Perform simple Monte Carlo using path discretization with \(m=50\) in order to approximate \[ A(0, T)=\frac{1}{T} \int_0^T S(t) d t \] by \(\tilde{A}(0, T)=\frac{1}{m} \sum_{j=1}^m S\left(t_i\right)\). Use \(n=10,000\) paths and report the standard error of the estimate.

Z=matrix(rnorm(n*m), n, m)
S=S0*exp(t(apply((r-v^2/2)*Dt+Z*v*sqrt(Dt),1,cumsum))) 
A=rowMeans(S)
payoff_5.a=exp(-r*M)*pmax(A-K,0)
(mean_5.a=mean(payoff_5.a))

## [1] 1.493

(se_5.a=sd(payoff_5.a)/sqrt(n))

## [1] 0.04219

(b) control variate

Repeat part (a) using the price of the geometric Asian call as a control variate, and report the new price estimate and standard error. For implementing the control variate, the exact price of the geometric Asian option is 1.271969.

EX=1.271969 # Control variate - exact Geom. Asian call price
A.geom=exp(rowMeans(log(S)))
X=exp(-r*M)*pmax(A.geom-K,0) # control variate values
Y=payoff_5.a
(rho_XY=cor(X,Y))

## [1] 0.9987

(b=sd(Y)/sd(X)*rho_XY)

## [1] 1.081

(mean_5.b=mean_5.a - b*(mean(X)-EX))

## [1] 1.405

(se_5.b=se_5.a*sqrt((1-rho_XY^2)))

## [1] 0.002139

(c) importance sampling

Repeat part (a) using importance sampling. Change the probability measure so that \(E\left[S_T\right]=K\), and experiment with at least two more measure changes of your choice. Report the new price estimates and standard errors.

MEANS=c(80,90,100)
mean_5.c=se_5.c=rep(0,length(MEANS))
for(i in 1:length(MEANS)){
  new.r=log(MEANS[i]/S0)/M
  ZZ=(new.r-v^2/2)*Dt+Z*v*sqrt(Dt)
  S=S0*exp(t(apply(ZZ,1,cumsum))) 
  A=rowMeans(S)
  log.ratio=dnorm(ZZ,(r-v^2/2)*Dt,v*sqrt(Dt),log=TRUE)-dnorm(ZZ,(new.r-v^2/2)*Dt,v*sqrt(Dt),log=TRUE)
  prod.ratio=exp( rowSums(log.ratio) )
  payoff=exp(-r*M)*pmax(A-K,0)*prod.ratio
  mean_5.c[i]=mean(payoff)
  se_5.c[i]=sd(payoff)/sqrt(n)
}

cbind(mean_5.c, se_5.c)

##      mean_5.c  se_5.c
## [1,]    1.480 0.02973
## [2,]    1.482 0.02138
## [3,]    1.494 0.01903

(a) [5 points]

Use the optimize () function to find the Kelly fraction \(f^*\) to invest, i.e. value of \(f\) which maximizes the expression in Q1.(g).

E_log_wealth = function(x){
  return( (2+1/x)*log(1+2*x) - (1/x-1)*log(1-x) - 1 ) 
}
f_ = optimize(E_log_wealth, c(0, 1), tol = 0.0001, maximum = T)$maximum
x = seq(0,1,.01)
plot(x, E_log_wealth(x), type ="l" )
points( f_, E_log_wealth(f_), pch = 16)

(b) [10 points]

Simulate 1,000 wealth paths consisting of 30 investments steps each using the Kelly fraction from the previous part (or, if you didn’t get it, just \(f=.7\) ). On the same axes, plot the mean simulated wealth, and the \(5 \%\) and \(95 \%\) quantiles of the simulated wealth at each step.

set.seed(123)
n=30; m=1000
R = matrix( runif(n*m, -1,2), n, m)
V = apply( 1+f_*R, 2, cumprod )

Vavg = apply( V, 1, mean )
Vlo = apply( V, 1, quantile, prob = .05 )
Vhi = apply( V, 1, quantile, prob = .95 )
plot( Vavg, type = "l", ylim = c(.001,25000) )
lines( Vlo, col=2 ); lines( Vhi, col = 2 )

# in log-scale
plot( Vavg, type = "l", log = "y", ylim = c(.001,25000) )

(c) [10 points]

In class, we mentioned that the Kelly criterion strategy reaches any wealth level faster than other strategies. You will test this claim with a simulation experiment. Fix a target wealth of 10,000 and simulate 1,000 paths with as many steps necessary to reach or exceed this level, under both the Kelly strategy \(\left(f^*\right)\) and the all-in strategy \((f=1)\). Report the average number of steps required to reach or exceed the target wealth under each strategy. Does it confirm the claim?

T_Kelly = T_all_in = integer(1000)

set.seed(123)
for( i in 1:1000 ){
  # for Kelly 
  V = 1; j = 0
  while( V < 10000 ){
    V = V * (1 + f_ * runif(1,-1,2) )
    j = j + 1
  }
  T_Kelly[i] = j
  
  # for all-in
  V = 1; j = 0
  while( V < 10000 ){
    V = V * (1 + runif(1,-1,2) )
    j = j + 1
  }
  T_all_in[i] = j
}

mean(T_Kelly)

## [1] 54.31

mean(T_all_in)

## [1] 99.17

The Kelly criterion takes a shorter time, on average, to reach the goal

(a) [5 points]

Create a scatter-plot of the values of the two industry factors. Do the factors look independent? Explain your answer based on what you know about asset behaviour.

plot( f[,1], f[,2] )

The two factors are positively correlated. This is expected, since we know that all asset returns have positive correlation, through their dependence on the overall “market” behavior.

(b) [3 points]

Use the sample covariance of the returns ( \(\mathrm{R}\) ) to calculate the minimum-variance portfolio weights.

S = var(R)
Si = solve(S)
w_mv_sample = rowSums(Si)/sum(Si)

(c) [10 points]

Use the industry factor model to estimate the return covariance matrix (from \(\boldsymbol{\beta} \boldsymbol{\Sigma}_F \boldsymbol{\beta}^{\top}+\boldsymbol{S}_e\), and use that to calculate the minimum-variance-portfolio weights.

S_e = diag( apply( e, 2, var ) )
S_f = cov(f) 

S_model = S_e + b %*% S_f %*% t(b)
Si = solve( S_model )
w_mv_model = rowSums(Si)/sum(Si)

(d) [3 points]

Use the weights you found in the last two parts to create the returns of each portfolio, and calculate their Sharpe ratio assuming a risk-free rate of 0.

R_mv_sample = R %*% w_mv_sample
(Sharpe_mv_sample = mean( R_mv_sample ) / sd( R_mv_sample ))

## [1] 0.1061

R_mv_model = R %*% w_mv_model
(Sharpe_mv_model = mean( R_mv_model ) / sd( R_mv_model ))

## [1] 0.1321

(e) [4 points]

Plot the cumulative gross returns of the two portfolios on the same axes.

plot( cumprod( 1 + R_mv_model ), type ="l", col = 2 , 
      ylab = "gross return", xlab = "day")
lines( cumprod( 1 + R_mv_sample ), col = 1  )
legend( "topleft", col=c(1,2), lwd=2, legend = c("Sample covariance", "Industry factor model covariance"))

(a) [5 points]

Calculate the log-returns of the ETFs, and plot the price and return series for the first ETF (DVEM).

plot(S[[1]])

r = lapply(S, FUN=function(x){diff(log(x))}) # log return
R = lapply(r, FUN=function(x){exp(x)-1}) # return
plot(R[[1]])

(b) [5 points]

Use factanal() to fit a 2-factor model to the correlation matrix of the returns. Report the factor loadings and idiosyncratic variances of you model.

Rmat = simplify2array(R)
fmod = factanal( Rmat, factors = 2,  lower = 0.005)
b = fmod$loadings
v = fmod$uniquenesses

(c) (10 points]

Simulate 250 daily log-returns using a multivariate Normal distribution with parameters given by the sample means and variances of the ETFs, and correlation matrix given by the previous factor model. Calculate and plot the cumulative net-returns of an equally weighted portfolio over the 10 ETFs.

MU = sapply(R, mean)
SD = sapply(R, sd)
VC = ( b%*%t(b) + diag(v) ) * (SD %*% t(SD)) 

Rsim = mvtnorm::rmvnorm(250, MU, VC)

R_eqwt = rowMeans( exp( apply(Rsim,2,cumsum) ) - 1 )
plot(R_eqwt, type = "l");

(a) [5 points]

Use fExoticOptions: : SimpleChooserOption () to find the exact price of the option.

# library(fExoticOptions)
# SimpleChooserOption(S = S0, X = K, time1 = 1, Time2 = 2, 
#                     r = r, b = r, sigma = v)@price

(b) [10 points]

Perform a simulation with \(n=10,000\) paths for pricing the chooser option. For each path, generate two prices: \(S_{T_1}\) at time \(T_1\) for determining the form of the payoff, and \(S_{T_2}\) at time \(T_2\) for determining the value of the payoff. Report the estimated price and its standard deviation.

set.seed(123); 
S0 = K = 100; T1 = 1; T2 = 2; r = .05; v = .2
n = 10000; dT = T2-T1
Z1 = rnorm(10000); Z2 = rnorm(10000)
S1 = S0*exp( (r-v^2/2)*T1 + v*sqrt(T1)*Z1 )
S2 = S1*exp( (r-v^2/2)*dT + v*sqrt(dT)*Z2 )
payoff.b = exp(-r*T2) * ifelse( S1 > K * exp(-r*dT), 
                              pmax( S2 - K, 0),
                              pmax( K - S2, 0) )
mean(payoff.b) 

## [1] 19.97

sd(payoff.b)/sqrt(n)

## [1] 0.2165

(c) [5 points]

Repeat the simulation experiment, but this time use the the Black Scholes formula to value the chosen option at time \(T_1\). In other words, you don’t have to simulate \(S_{T_2}\); just simulate \(S_{T_1}\) and calculate the exact value of the chosen call or put option at time \(T_1\) using Black-Scholes. Report the estimated price and its standard deviation (should be smaller). (Hint: use f0ptions : : GBSOption () to get the Black-Scholes price.)

# library(fOptions)
# payoff.c = exp(-r*T1) * ifelse( S1 > K * exp(-r*dT), 
#                                 GBSOption( "c", S1, K, dT, r, r, v )@price,
#                                 GBSOption( "p", S1, K, dT, r, r, v )@price )
# mean(payoff.c) 
# sd(payoff.c)/sqrt(N)

(d) [10 points]

Repeat part (b), using the following dynamics for the asset: \[ d S_t=r S_t d t+\sigma \log \left(S_t\right) d W_t \] Note that the price process no longer follows GBM, so use path discretization with \(m=20\) steps per path to approximate \(S_{T_1}\) and \(S_{T_2}\). Report the estimated price and its standard deviation. (Note: the option choice rule at \(T_1\) does not change with asset price dynamics, since put-call parity holds for any model.)

m = 20; Dt = T2/m
Z = matrix( rnorm( n * m ), n, m)
S = matrix( S0, n, m + 1 )
# iterate over all paths
for(i in 1:m){
  S[,i+1] = S[,i]  + r * S[,i] * Dt + v * log(S[,i]) * sqrt(Dt) * Z[,i] 
}
payoff.d = exp(-r*T2) * ifelse( S[, m/2 + 1] > K * exp(-r*dT),
                                pmax( S[, m + 1] - K, 0),
                                pmax( K - S[, m + 1], 0) )
mean(payoff.d)

## [1] 9.519

sd(payoff.d) / sqrt(n)

## [1] 0.01262